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Stability as a natural selection mechanism on interacting networks 

Juan I. Perottij 11121 * Orlando V. BilloniJ 1113 Francisco A. Tamaritj 1113 

Sergio A. Cannas^^ 

Biological networks of interacting agents exhibit similar topological properties for a wide 
range of scales, from cellular to ecological levels, suggesting the existence of a common 
evolutionary origin. A general evolutionary mechanism based on global stability has been 
proposed recently [J I Perotti, et al., Phys. Rev. Lett. 103, 108701 (2009)]. This mech- 
anism was incorporated into a model of a growing network of interacting agents in which 
each new agent's membership in the network is determined by the agent's effect on the 
network's global stability. In this work, we analyze different quantities that characterize 
the topology of the emerging networks, such as global connectivity, clustering and average 
nearest neighbors degree, showing that they reproduce scaling behaviors frequently ob- 
served in several biological systems. The influence of the stability selection mechanism on 
the dynamics associated to the resulting network, as well as the interplay between some 
topological and functional features are also analyzed. 



I. Introduction 

The concept of networks of interacting agents has 
proven, in the last decade, to be a powerful tool 
in the analysis of complex systems (for reviews, 
see Refs. [1-4]). Although not new, with the ad- 
vent of high performance computing, this theoreti- 
cal construction opened a new door for the statisti- 
cal physics methodology in the analysis of systems 
composed by a large number of units that interact 
in a complicated way. This allowed to get new in- 
sights about the dynamical behavior of systems as 
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complex as biological and social systems. In ad- 
dition, it constitutes a basic backbone upon which 
relatively simple models can be constructed in a 
bottom-up strategy. 

As a modeling tool, the definition of an interac- 
tion network for a given system is frequently not 
unique (see for example the case of protein-protein 
interaction networks [5-7]), depending on the coarse 
grain level of the approach. Nevertheless, many 
topological properties appear to be independent of 
the definition of the network. Moreover, some of 
those properties have emerged in the last years as 
universal features among systems otherwise consid- 
ered very different from each other. In particular, 
the following properties are characteristic of most 
biological networks, (a) Small worldness: all of 
them exhibit high clustering Cc and relatively short 
path length L, compared with random networks. L 
is defined as the minimum number of links needed 
to connect any pair of nodes in the network and Cc 
is defined as the fraction of connections between 
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topological neighbors of any sitepQ. (b) Scale free 
degree distribution: the degree distribution P{k) 
(the probability of a node to be connected to k 
other ones) presents a broad tail for large values 
of k. In some cases, the tail can be approached 
by a power law P(k) ~ fc~ 7 with degree exponents 
7 < 3 for a wide range of scales, while in others, a 
cutoff appears for some maximum degree k max ] in 
the latter, the degree distribution is generally well 
described by P(k) - fe-7 e - fe A— [1, 7-12]. In 
any case, the networks present a nonhomogeneous 
structure, very different from that expected in a 
random network, (c) Scaling of the clustering co- 
efficient: in many natural networks, it is observed 
that the clustering coefficient of a node with degree 
k follows the scaling law Cc(k) ~ k~@ , with (3 tak- 
ing values close to one. This has been interpreted 
as an evidence for a modular structure organized in 
a hierarchical way [T3]. (d) Disassortative mixing 
by degree: in most biological networks, highly con- 
nected nodes tend to be preferentially connected to 
nodes with low degree and vice versa [Tl] . 

These properties are observed for a wide range 
of scales, from the microscopic level of genetic, 
metabolic and proteins networks to the macro- 
scopic level of communities of living beings (eco- 
logical networks). Such ubiquity suggests the ex- 
istence of some natural selection process that pro- 
motes the development of those particular struc- 
tures [2j. One possible constraint general enough 
to act across such a range of scales is the proper 
stability of the underlying dynamics. 

Growing biological networks involve the coupling 
of at least two dynamical processes. The first one 
concerns the addition of new nodes, attached dur- 
ing a slow evolutionary (i.e., species lifetime) pro- 
cess. A second one is the node dynamics which 
affects and in turn is affected by the growing pro- 
cesses. It is reasonable to expect that the network 
topologies we finally witness could have emerged 
out of these coupled processes. Consider, for ex- 
ample, the case of an ecological network like a food 
web, where nodes are species within an ecosystem 
and edges are consumer-resource relationships be- 
tween them. New nodes are added during evo- 
lutionary time scales, through speciation or mi- 
gration of new species. Then, the network grows 
through community assembly rules, strongly influ- 
enced by the underlying dynamics of species and 
specific interactions among them [151 116] . The con- 



sequence of adding a new member with a given 
connectivity affecting a global in/stability, is repre- 
sented in this case by the aboundance/lack of food 
[T7] . Notice that each new member may not only 
result in its own addition/rejection to the system, 
but it can also promote avalanches of extinctions 
amongst existing members. 

The above ingredients were recently incorporated 
into a simple model of growing networks under sta- 
bility constraints [H]. Numerical simulations on 
this model showed that, indeed, complex topology 
can emerge out of a stability selection pressure. 
In the present work, we further explore different 
topological and dynamical properties predicted by 
the model, whose definition is reviewed in section 
TT7I The results are presented in sections IIII.I and 
IV. I In section IIII.I we analyze the topological fea- 
tures that emerge in growing networks under sta- 
bility constraint. In section ITVl we show that this 
constraint not only induces topological features of 
the resulting networks but also influences the as- 
sociated dynamics. A discussion of the results is 
presented in section IVTI 



II. The Model 

Let us consider a system of n interactive agents, 
whose dynamics is given by a set of differen- 
tial equations dx/dt = F(x), where x is an n- 
component vector describing the relevant state 
variables of each agent and F is an arbitrary non- 
linear function. One could imagine that x in differ- 
ent systems may represent concentrations of some 
hormones, the average density populations in a 
food web, the concentration of chemicals in a bio- 
chemical network, or the activity of genes in a gene 
regulation net, etc. We assume that a given agent 
i interacts only with a limited set of fe, < n other 
agents; thus, Fi depends only on the variables be- 
longing to that set. This defines the interaction 
network. 

We assume that there are two time scales in the 
dynamics. Let f m be the average frequency of the 
incoming flux of new agents (migration, mutation, 
etc.). This defines a characteristic time r m = /^ • 
On the long time scale t 3> r m (much larger than 
the observation time) new agents arrive to the sys- 
tem and start to interact with some of the previous 
ones. Some of them can be incorporated into the 
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system or not, so n (and the whole set of differential 
equations) can change. Once a new agent starts to 
interact with the system, we will assume that the 
enlarged system evolves towards some stationary 
state with characteristic relaxation time r re i <C r m • 
Then, in the short time scale r re i < f « T m we 
can assume that n is constant and the dynamics 
already led the system to a particular stable sta- 
tionary state x* defined by F(x*) — 0. Following 
May's ideas [18] , we assume that the only attractors 
of the dynamics are fixed points. Nevertheless, the 
proposed mechanism is expected to work, as well, 
for more complex attractors (e.g, limit cycles). 

The stability of the solution x* is determined by 
the eigenvalue with maximum real part of the Ja- 
cobian matrix 



**ij 



m 

dxj 



(1) 



A new agent will be incorporated to the network 
if its inclusion results in a new stable fixed point, 
that is, if the values of the interaction matrix Ojj 
are such that the eigenvalue with maximum real 
part A m of the enlarged Jacobian matrix is nega- 
tive (A m < 0). Assuming that isolated agents will 
reach stable states by themselves after certain char- 
acteristic relaxation time, the diagonal elements of 
the matrix a^ are negative and given unity value 
to further simplify the treatment [IS] ■ The inter- 
action values, (i.e., the non-diagonal matrix ele- 
ments ctij) will take random values (both positive 
and negative) taken from some statistical distri- 
bution. In this way, we have an unbounded en- 
semble of systems [TB] characterized by a "growing 
through stability" history. Randomness would be 
self-generated through the addition of new agents 
processes. Each specific set of matrix elements, af- 
ter addition, defines a particular dynamical system 
and the subsequent analysis for time scales between 
successive migrations is purely deterministic. 

The model is then defined by the following algo- 
rithm 19 . At every step, the network can either 
grow or shrink. In each step, an attempt is made 
to add a new node to the existing network, starting 
from a single agent (n = 1). Based on the stability 
criteria already discussed, the attempt can be suc- 
cessful or not. If successful, the agent is accepted, 
so the existing n x n matrix grows its size by one 
column and one row. Otherwise, the novate agent 



will have a probability to be deleted together with 
some other nodes, as further explained below. 

More specifically, suppose that we have an al- 
ready created network with n nodes, such that the 
n x n associated interaction matrix ajj is stable. 
Then, for the attachment of the (n + l) t h node, 
we first choose its degree k n +i randomly between 
1 and n with equal probability. Then, the new 
agent interaction with the existing network member 
i is chosen, such that non-diagonal matrix elements 
(Oi j7l 4-i, On+i.i) (i = 1, . . . , n) are zero with proba- 
bility 1 — fc n +i/n and different from zero with prob- 
ability fc n _|_i/n; to each non-zero matrix element 
we assign a different real random value uniformly 
distributed in [—&,&]. b determines the interaction 
range variability and it is one of the two parameters 
of the model [20] . 

Then, we calculate numerically A m for the result- 
ing (n + 1) x (n + 1) matrix. If A m < 0, the new 
node is accepted. If A m > 0, it means that the 
introduction of the new node destabilized the en- 
tire system and we will impose that, the new agent 
is either eliminated or it remains but produces the 
extinction of a certain number of previous existing 
agents. In order to further simplify the numerical 
treatment, we allow up to q < k n +\ extinctions, 
taken from the set of fc n +i nodes connected to the 
new one; q is the other parameter of the model. To 
choose which nodes are to be eliminated, we first 
select one with equal probability in the set of k n +i 
and remove it. If the resulting n x n matrix is sta- 
ble, we start a new trial; otherwise, another node 
among the remaining fc n +i — 1 is chosen and re- 
moved, repeating the previous procedure. If after 
q removals the matrix remains unstable, the new 
node is removed (we return to the original n x n 
matrix and start a new trial). The process is re- 
peated until the network reaches a maximum size 
n = n max (typically n max = 200) and restarted M 
times from n = 1 to obtain statistics of the net- 
works (typically M = 10 5 ). 

III. Topological properties 

i. Connectivity 

First, we analyzed the average connectivity C(n), 
defined as the fraction of non-diagonal matrix el- 
ements different from zero, averaged over different 
runs. In Fig. [TJ we show the typical behavior of 
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Figure 1 : Connectivity as a function of the network 
size for q — 3, n max = 200 and different values of b. 
The symbols correspond to numerical simulations 
and the dashed lines to power law fittings of the 
tails C(n) = B n~( 1+e \ The insets show the fitting 
values B and e as a function of b 



C(n) for different values of b (we found that C(n) 
is completely independent of q). The connectivity 
presents a power law tail for large values of n. From 
a fitting of the tail with a power law (see insets in 
Fig. [1} we obtain the scaling behavior 



C(n) 



i+<0 



(2) 



for large values of n, where a is the variance of the 
non-diagonal elements of the stability matrix (a — 
b 2 /3 for the uniform distribution) and uj — 0.7±0.1. 
From the inset of Fig. [TJ we see that the exponent 
e shows a weak dependency on 6, taking values in 
the range (0.1,0.3) . It is interesting to compare 
Eq. ((2J with May's stability line for random net- 
works pj| C(n) — (an)^ 1 . It is easy to see that 
Eq. © lies above May's stability line for network 
sizes up to ~ 10 6 [3T]. This shows that networks 
growing under stability constraint develop partic- 
ular structures whose probability in a completely 
random ensemble is almost zero. In other words, 
the associated matrices belong to a subset of the 
random ensemble with zero measure and therefore 
they are only attainable through a constrained de- 
velopment process. In the next sections, we explore 
the characteristics of those networks. 

In Fig. [3J we plotted the connectivity for differ- 
ent biological networks across three orders of mag- 
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Figure 2: Connectivity as a function of the network 
size for different biological networks. The straight 
line is a power law fitting C(n) = An^ l - 1+e \ giv- 
ing an exponent e = 0.2 ± 0.1 (R 2 = 0.92). Data 
extracted from: [5[ [7] (protein-protein interaction 
networks); [5] (metabolic networks); [2H [53] (food 
webs). 



nitude of network size scales, using data collected 
from the literature. We see that the data are very 
well fitted by a single power law C(n) ~ n , in 
a nice agreement with the average value e = 0.2 
predicted by the present model. It is worth men- 
tioning that the behavior C(n) ~ n~( 1+e ' has also 
been obtained in a self organized criticality model 
of Food Webs [H]. 

ii. Degree distribution 

The degree distribution P(k) of the network was 
analyzed in detail in Ref. [JJj]. We briefly summa- 
rize the main results here. In Fig. [3j we illustrate 
the typical behavior of P(k). It presents a power 
law tail P(k) ~ fc" 7 for values of k > 20, with a 
finite size drop at k = n max . The degree exponent 
7 takes values between 2 and 3 for values of b in 
the interval 6e (1.5,3.5), which become almost in- 
dependent of q as it increases. The exponent 7 can 
also fall below 2 when the global stability constraint 
is replaced by a local one. The qualitative struc- 
ture of P(k) remains when the stability criterium 
A m < is relaxed by the condition X m < A, with 
A some small positive number. In other words, 
the power law tail emerges also when the addition 
of new nodes destabilizes the dynamics, provided 
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Figure 3: Degree distribution P(k) for q = 3, 
n ma x — 200 and different values of 6; the dashed 
lines correspond to power law fittings of the tail 
P{k) ~ fc~ 7 . Logarithmic binning has been used 
to smooth the curves. 
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Figure 4: Average network size as a function of the 
time measured in number of trials for 6 = 2 and 
q = 3. The continuous red line corresponds to the 
power law i 1 / 2 , 



that the characteristic time to leave the fixed point 
t = A^ 1 is large enough to become comparable to 
the migration time scale r m [19] . 



iii. Network growth and clustering proper- 
ties 

Networks grown under stability constraint also dis- 
play small world properties. The average clustering 
coefficient decays with the network size as Cc(n) ~ 
(which is slower than the 1/n decay in a 



-0.75 



random net), while the average path length L be- 
tween two nodes increases as L(n) ~ A In (n + C) 
[T9"] . A similar behavior is observed in the Barabasi- 
Albert model [1], where the clustering can be ap- 
proximated by a power law with the same exponent, 
although the exact scaling is [27] Cc(n) ~ (In n) 2 /n 
(therefore that behavior cannot be excluded in the 
present model). While this suggests the presence of 
an underlying preferential attachment rule mecha- 
nism, a detailed analysis has shown that this is not 
the dominant mechanism |19j . The behavior of Cc 
and L is linked with the selection dynamics ruling 
which node is accepted or rejected. The stability 
constraint favors the nodes with few links, since 
they modify the matrix a^j stability much less than 
new nodes with many links (of course this is re- 
flected in the P(k) density). Thus, most frequently, 



the network grows at the expense of adding nodes 
with one or few links, producing an increase of L 
and a decrease of Cc, but sporadically, a highly con- 
nected node is accepted, decreasing L and increas- 
ing Cc(n) [15]. Those fluctuations lead to a slow 
diffusive-like growth of the network size n(t) ~ t 1 ^ 2 
(See Fig. |H). 

Another quantity of interest is the average clus- 
tering Cc(k) as a function of the degree k. A typical 
example is shown in Fig. [SJ We see that Cc{k) de- 
creases monotonously with k and displays a power 
law tail Cc(k) ~ k~@ with an exponent (3 ss 0.9, 
close to one. The exponent appears to be com- 
pletely independent of b and q. This behavior is 
indicative of a modular structure with hierarchical 
organization [13]. Notice that this power law de- 
cay appears for degrees k > 20, precisely the same 
range of values for which the degree distribution 
P(k) displays a power law tail (see subsection HT7|) . 

iv. Mixing by degree patterns 

To analyze the mixing by degree properties of the 
networks selected by the stability constraint, we 
calculated the average degree k nn among the near- 
est neighbors of a node with degree k. In Fig. \§\ we 
see that k nn decays with a power law k nn ~ k 
for k > 20, with an exponent 6 close to —0.25, in 
a clear disassortative behavior. This result is also 
consistent with previous works showing that assor- 
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Figure 5: Average clustering coefficient Cc(k) as a 
function of the degree k for b = 2, q = 3 and dif- 
ferent values of n max . The dashed black line corre- 
sponds to the power law fc~ 09 . 
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Figure 6: Average nearest neighbors degree k nn (k) 
of a node with degree k for q = 3 and different 
values of b and n ma£C . The dashed line corresponds 
to the power law fc~ 0,25 . 



tative mixing by degree decreases the stability of 
a network, i.e., the maximum real part A m of the 
eigenvalues of random matrices of the type here 
considered increases faster on assortative networks 
than on disassortative ones [29] . 



IV. Dynamical properties 

In the previous section, we analyzed different topo- 
logical properties that are selected by the stabil- 
ity constraint, i.e., properties associated to the un- 
derlying adjacency matrix, regardless of the values 
of the interaction strengths. We now analyze the 
characteristics of the dynamics associated to the 
networks emerging from such constraint. In other 
words, we investigate the statistics of values of the 
non null elements a%j ^ 0. 

First of all, we calculate the probability distribu- 
tion of values for a single non null matrix element 
Ojj of the final network with size n = n max . The 
typical behavior is shown in Fig. [7) We see that 
Piflij) is an even function, almost uniform in the 
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Figure 7: Probability density of the matrix ele- 



ments a™ for b = 2, q = 3 and n T< 
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tively enhances the stability of random matrices 
|31j . In an ecological network, this typically cor- 
responds to a predator-prey or parasite-host inter- 
action. To check for the presence of such type of 
interval [-b,b], with a small cusp around ay = 0. interactions, we calculated the correlation {a tj a jt ), 



This shows that stability is not enhanced by a par- 
ticular sign or absolute value of the individual in- 
teraction coefficients. It has been shown recently 
that the presence of anticorrclated links between 
pairs of nodes (i.e., links between pairs of nodes 
(i,j) such that sign(aij) — —sign(aji)) significa- 



where the average is taken over pairs of nodes with 
a double link (ay ^ and dji ^ 0). 

In Fig. [51 we show (o»jOj») as a function of the 
network size n. We see that this correlation is neg- 
ative for any value of n and saturates into a value 
{dijdji) ?s —0.65 for large values of n. In the in- 
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Figure 8: Correlation lunction (aijCLji) between a 
pair of nodes (i,j) with a double link as a function 
of the network size, for b = 2 and q = 3. The aver- 
age was calculated over all pair of sites with dou- 
ble link in networks with the same size. The inset 
shows a comparison between the fraction of double 
links in the present network (ry) (green circles) and 
a random network of the same size and connectivity 
C(n): (j])ran = C / (2 — C) (continuous line). Yel- 
low triangles correspond to a numerical calculation 
of (rj)ran- The dashed line corresponds to a power 
law n" - 68 . 



set of Fig. [SI we compare the average fraction of 
double links (77) with the corresponding quantity 
for a completely random network with the same 
connectivity C(n), that is, a network where all 
edges are independently distributed with probabil- 
ity P(ajj 7^ 0) = C(n). Then, the probability of 
having a link between an arbitrary pair of sites is 
Pd = C(n)(2 — C{n)) and the average degree per 
node (fc) = pan. Hence 



yljran — 



C 2 (n)n _ C(n) 
(k) ~ 2 - C{n) 



(3) 



Then for large values of n, we have (f]} ran ~ C(n) ~ 
n -(i+e)_ F rom the inset of Fig. we see that 
(rf) ~ n~ - 68 when n > 1 in the present case. The 
fraction of double links is considerable larger than 
in a random network. The two results of Fig. [^to- 
gether show that the present networks have indeed 
a significantly large number of anticorrelated pair 
interactions. 

Next, we calculated the correlation 
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Figure 9: Correlation function {a-ijCLji) / {\aij\) 2 be- 
tween the matrix elements linking a node i and its 
neighbors j as a function of its degree fcj, for b = 2, 
q = 3 and different values of n max . The inset shows 
the average fraction of anticorrelated links (k) as a 
function of ki. 



(aijdji) / (\aij\) 2 between the matrix elements, 
linking a node i and its neighbors j, as a function 
of its degree ki, where the average is taken only 
on the double links. From Fig. [HI we see that 
the absolute value of the correlation presents a 
maximum around ki = 25 and tends to zero as the 
degree increases. The inset of Fig.[9]shows that the 
average fraction of anticorrelated links (k) (i.e., 
# anticorrelated links/total # double links) tends 
to 1/2 as the degree increases. We can conclude 
from these results that the interactions strengths 
between the hubs and their neighbors are almost 
uncorrelated. This suggests that the influence of 
hubs in the stabilization of the dynamics is mainly 
associated to their topological role (e.g., reduction 
of the average length L) rather than to the nature 
of their associated interactions. 



V. Discussion 

The recent advances in the research on networks 
theory in biological systems have called for a deeper 
understanding about the relationship between net- 
work structure and function, based on evolutionary 
grounds [3]. In this work, we have shown that a 
key factor to explain the emergency of many of the 
complex topological features commonly observed in 
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biological networks could be just the stability of the 
underlying dynamics. Stability can then be consid- 
ered as an effective fitness acting in all biological 
situations. The results presented in Fig. [2] for the 
connectivity of real biological networks at different 
network size scales support this conclusion. In ad- 
dition, the present approach (although based on a 
very simple model) allows to draw some conclusions 
about the interplay between network structure and 
function that could be of general applicability. The 
present results suggest that hubs play mainly a 
topological role of linking modules (disassortativ- 
ity, low clustering, uncorrelated links), while low 
connected nodes inside modules enhance stability 
through the presence of many anticorrelated inter- 
actions. 

The stabilizing effects of some of the topological 
and functional network features here analyzed have 
been previously addressed separately (small world 
[52] . dissasortative mixing [2j|] EU], anticorrelated 
interactions [31]). However, the present analysis 
suggests that the simultaneous observance of all of 
them is highly unlikely to be a result of a purely 
random process. Such delicate balance of specific 
topological and functional features would only be 
attainable through a slow, evolutionary stability se- 
lection process. 

In particular, the above scenario agrees very well 
with the observed structures in cellular networks. 
For instance, the scaling behavior of Cc(k), dis- 
played in Fig. [5] has been observed in metabolic 
[28] and protein [6] [7] networks. Disassortative 
mixing by degree is another ubiquous property of 
those systems and indeed a very similar behavior 
to that shown in Fig. [5] has been observed in cer- 
tain protein-protein interaction networks [7. Also, 
the available data for the degree distribution in all 
those cases are consistent with a power law behav- 
ior with an exponent 7 between 2 and 2.5 [Tj. The 
agreement with the whole set of properties pre- 
dicted by the model suggests that stability could 
be a key evolutionary factor in the development of 
cellular networks. 

The situation is a bit different in the case of eco- 
logical networks, where the predictions of the model 
do not completely agree with the observations, spe- 
cially those related to food webs. On the one hand, 
food webs usually display also disassortative mix- 
ing by degree, modularity and relatively low small 
worldness 33 (rather low values of clustering, com- 



pared with other biological networks) , in agreement 
with the present predictions. Regarding the scal- 
ing behavior of C(n) [23], this is the topic of an old 
debate in ecology (see Dunne's review in Ref. [23] 
for a summary of the debate). While in general 
it is expected a power law behavior, the value of 
the exponent (and the associated interpretations) 
is controversial, due to the large dispersion of the 
available data, the rather small range of network 
sizes available and, in some cases, the low resolu- 
tion of the data [25] . The consistency of the scal- 
ing shown in Fig. [2] for a broad range of size scales 
suggests that the ecological debate should be recon- 
sidered in a broader context of evolutionary growth 
under stability constraints. 

On the other hand, the degree distribution of 
food webs is not always a power law, but it fre- 
quently exhibits an exponential cutoff at some max- 
imum characteristic degree k max [9J [TU] ■ Such vari- 
ance between food webs and other biological net- 
works is probably related to the way ecosystems 
assemble and evolve compared with other systems. 
While the hypothesis of the present model are gen- 
eral enough to apply in principle to any biological 
system, that difference suggests that stability is not 
enough to explain the observed structures in food 
webs, but further constraints should be included 
to account for them. For instance, at least two 
different (although closely related) constraints are 
known that can generate a cutoff in the degree dis- 
tribution: aging and limited capacity of the nodes 
[34] . In the former, nodes can become inactive with 
some probability through time (in the sense that 
they stop interacting with new agents), while in 
the latter they systematically pay a "cost" every 
time a new link is established with them, so that 
they become inactive when some maximum degree 
is reached. One can easily imagine different situ- 
ations in which mechanisms of that type become 
important in the evolution of ecological webs, ei- 
ther by limitations in the available resources or by 
dynamical changes in the diet of species due to ex- 
ternal perturbations (for instance, there are many 
factors that constrain a predator's diet; see Ref. 
[S] and references therein for a related discussion). 
Mechanisms of these kind can be easily incorpo- 
rated into the model, serving as a base for the de- 
scription of more complex behaviors in particular 
systems like food webs. 

Finally, it would be interesting to analyze the 
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relationship between dynamical stability in evolv- 
ing complex networks and synchronization, a topic 
about which closely related results have been re- 
cently published [55] . 
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